// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//                     avtMatrixMarketFileFormat.C                           //
// ************************************************************************* //

#include <avtMatrixMarketFileFormat.h>

#include <string>

#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkUnstructuredGrid.h>

#include <avtDatabaseMetaData.h>
#include <avtDatabase.h>

#include <Expression.h>

#include <InvalidVariableException.h>
#include <InvalidFilesException.h>

extern "C" {
#include "mmio.h"
}

using     std::string;


// ****************************************************************************
//  Method: avtMatrixMarketFileFormat constructor
//
//  Programmer: Jeremy Meredith
//  Creation:   March 13, 2009
//
// ****************************************************************************

avtMatrixMarketFileFormat::avtMatrixMarketFileFormat(const char *fn)
    : avtSTSDFileFormat(fn)
{
    filename = fn;
    dataread = false;
    width = height = 0;
    matrix = NULL;
}


// ****************************************************************************
//  Method: avtMatrixMarketFileFormat::FreeUpResources
//
//  Purpose:
//      When VisIt is done focusing on a particular timestep, it asks that
//      timestep to free up any resources (memory, file descriptors) that
//      it has associated with it.  This method is the mechanism for doing
//      that.
//
//  Programmer: Jeremy Meredith
//  Creation:   March 13, 2009
//
// ****************************************************************************

void
avtMatrixMarketFileFormat::FreeUpResources(void)
{
    if (matrix)
        matrix->Delete();
    width = height = 0;
    dataread = false;
}


// ****************************************************************************
//  Method: avtMatrixMarketFileFormat::PopulateDatabaseMetaData
//
//  Purpose:
//      This database meta-data object is like a table of contents for the
//      file.  By populating it, you are telling the rest of VisIt what
//      information it can request from you.
//
//  Programmer: Jeremy Meredith
//  Creation:   March 13, 2009
//
// ****************************************************************************

void
avtMatrixMarketFileFormat::PopulateDatabaseMetaData(avtDatabaseMetaData *md)
{
    ReadData();

    AddMeshToMetaData(md, "cellmesh", AVT_RECTILINEAR_MESH, NULL, 1,0, 2,2);
    AddScalarVarToMetaData(md, "cellvals", "cellmesh", AVT_ZONECENT);

    AddMeshToMetaData(md, "nodemesh", AVT_RECTILINEAR_MESH, NULL, 1,0, 2,2);
    AddScalarVarToMetaData(md, "nodevals", "nodemesh", AVT_NODECENT);
}


// ****************************************************************************
//  Method: avtMatrixMarketFileFormat::GetMesh
//
//  Purpose:
//      Gets the mesh associated with this file.  The mesh is returned as a
//      derived type of vtkDataSet (ie vtkRectilinearGrid, vtkStructuredGrid,
//      vtkUnstructuredGrid, etc).
//
//  Arguments:
//      meshname    The name of the mesh of interest.  This can be ignored if
//                  there is only one mesh.
//
//  Programmer: Jeremy Meredith
//  Creation:   March 13, 2009
//
// ****************************************************************************

vtkDataSet *
avtMatrixMarketFileFormat::GetMesh(const char *meshname)
{
    ReadData();

    int w=width;
    int h=height;

    // If we want the mesh as a cell-centered variable, we'll
    // just create the mesh one larger in each dimension
    if (string(meshname)=="cellmesh")
    {
        w++;
        h++;
    }

    // rectilinear grid
    vtkRectilinearGrid   *rgrid   = vtkRectilinearGrid::New(); 

    vtkFloatArray   *coords[3];
    int dims[3] = {w, h, 1};
    for (int i = 0 ; i < 3 ; i++)
    {
        coords[i] = vtkFloatArray::New();
        coords[i]->SetNumberOfTuples(dims[i]);

        for (int j = 0 ; j < dims[i] ; j++)
        {
            coords[i]->SetComponent(j, 0, j);
        }
    }
    rgrid->SetDimensions(dims);
    rgrid->SetXCoordinates(coords[0]);
    rgrid->SetYCoordinates(coords[1]);
    rgrid->SetZCoordinates(coords[2]);
    coords[0]->Delete();
    coords[1]->Delete();
    coords[2]->Delete();

    return rgrid;        

}


// ****************************************************************************
//  Method: avtMatrixMarketFileFormat::GetVar
//
//  Purpose:
//      Gets a scalar variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      varname    The name of the variable requested.
//
//  Programmer: Jeremy Meredith
//  Creation:   March 13, 2009
//
// ****************************************************************************

vtkDataArray *
avtMatrixMarketFileFormat::GetVar(const char *varname)
{
    ReadData();

    matrix->Register(NULL);
    return matrix;
}


// ****************************************************************************
//  Method: avtMatrixMarketFileFormat::GetVectorVar
//
//  Purpose:
//      Gets a vector variable associated with this file.  Although VTK has
//      support for many different types, the best bet is vtkFloatArray, since
//      that is supported everywhere through VisIt.
//
//  Arguments:
//      varname    The name of the variable requested.
//
//  Programmer: Jeremy Meredith
//  Creation:   March 13, 2009
//
//  Modifications:
//
//   Mark C. Miller, Wed Jun 17 17:46:18 PDT 2009
//   Replaced CATCHALL(...) with CATCHALL
//
//   Jeremy Meredith, Thu Jun 13 15:13:02 EDT 2013
//   Added symmetric support.  For dense arrays, the file only contains
//   values at or below the diagonal.  For sparse, we just write to 
//   both i,j and j,i.
//
// ****************************************************************************

vtkDataArray *
avtMatrixMarketFileFormat::GetVectorVar(const char *varname)
{
    return NULL;
}

void
avtMatrixMarketFileFormat::ReadData()
{
    if (dataread)
        return;

    FILE *f;
    if ((f=fopen(filename.c_str(), "r")) == NULL)
        EXCEPTION2(InvalidFilesException, filename, "File didn't exist.");

    TRY
    {
        // Get the matrix type
        MM_typecode matcode = "\n\n";
        if (mm_read_banner(f, &matcode) != 0)
            EXCEPTION2(InvalidFilesException, filename,
                       "Header was malformed.");

        // If we're on the mdserver, we're done and now assume a good
        // header is suffient
        if (!avtDatabase::OnlyServeUpMetaData())
        {
            // Read the types we support
            if (mm_is_matrix(matcode) &&
                (mm_is_real(matcode) ||
                 mm_is_integer(matcode)) &&
                mm_is_array(matcode) &&
                (mm_is_symmetric(matcode) ||
                 mm_is_general(matcode)))
            {
                bool symm = mm_is_symmetric(matcode);
                // Dense, real-valued, general/symmetric array
                mm_read_mtx_array_size(f, &height, &width);
                int size = width*height;
                matrix = vtkDoubleArray::New();
                matrix->SetNumberOfComponents(1);
                matrix->SetNumberOfTuples(size);
                double val;
                for (int j=0; j<width; j++)
                {
                    for (int i=0; i<height; i++)
                    {
                        if (symm && i<j)
                        {
                            val = matrix->GetComponent(j*width+i, 0);
                        }
                        else
                        {
                            int nread = fscanf(f, "%lg\n", &val);
                            if (nread != 1)
                                EXCEPTION2(InvalidFilesException, filename,
                                           "Premature EOF.");
                        }
                        matrix->SetComponent(i*width+j, 0, val);
                    }
                }
            }
            else if (mm_is_matrix(matcode) &&
                     (mm_is_real(matcode) ||
                      mm_is_integer(matcode)) &&
                     mm_is_coordinate(matcode) &&
                     (mm_is_symmetric(matcode) ||
                      mm_is_general(matcode)))
            {
                bool symm = mm_is_symmetric(matcode);
                // Sparse, real-valued, general/symmetric array
                int nvals;
                mm_read_mtx_crd_size(f, &height, &width, &nvals);
                int size = width*height;
                matrix = vtkDoubleArray::New();
                matrix->SetNumberOfComponents(1);
                matrix->SetNumberOfTuples(size);
                for (int e=0; e<size; e++)
                    matrix->SetComponent(e, 0, 0.0);
                for (int v=0; v<nvals; v++)
                {
                    int i, j;
                    double val;
                    int nread = fscanf(f, "%d %d %lg\n", &i, &j, &val);
                    --i; // 1-origin in the file, so
                    --j; // convert it to 0-origin
                    if (nread != 3)
                        EXCEPTION2(InvalidFilesException, filename,
                                   "Premature EOF.");
                    matrix->SetComponent(width*i+j, 0, val);
                    if (symm)
                        matrix->SetComponent(width*j+i, 0, val);
                }
            }
            else
            {
                // Not supported yet....
                EXCEPTION2(InvalidFilesException, filename,
                           "Unsupported matrix type.");
            }
        }
    }
    CATCHALL
    {
        fclose(f);
        RETHROW;
    }
    ENDTRY;

    fclose(f);
}
